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We study the Ruderman-Kittle-Kasuya-Yosida (RKKY) interaction in the presence of spin polar- 
ized two dimensional Dirac fermions. We show that a spin polarization along the z-axis mediates 
an anisotropic interaction which corresponds to a XXZ model interaction between two magnetic 
moments. For undoped graphene, while the x part of interaction keeps its constant ferromagnetic 
sign, its 2 part oscillates with the distance of magnetic impurities, R. A finite doping causes that 
both parts of the interaction oscillate with R. We explore a beating pattern of oscillations of the 
RKKY interaction along armchair and zigzag lattice directions, which occurs for some certain values 
of the chemical potential. The two characteristic periods of the beating are determined by inverse 
of the difference and the sum of the chemical potential and the spin polarization. 
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I. INTRODUCTION 

The charge and spin oscillatory interactions in metals 
has attracted considerable attention both on the theo- 
retical and experimental sides^l. Ruderman and Kittle^ 
suggested that the spin oscillatory interaction in met- 
als could provide a long-range interaction between nu- 
clear spins in metals. Afterwards, Kasuya and YosidaP 
extended the theory to include the long-range interac- 
tion between magnetic impurities and thus the combined 
refers to RKKY interaction. 

The recent discovery of graphene^, the two- 
dimensional crystal of carbon atoms, has provided a new 
material with a peculiar structure for the charge and 
the spin interactions. This stable crystal has already 
attracted considerable attention because of its unusual 
effective many-body properties^!] that follow from chi- 
ral nature of linearly dispersing low energy excitations 
described by pair of Dirac cones at the K and K' edges 
of the first Brillouin zone. 

The RKKY interaction in pristine graphene has been 
studied by several groups^^. Due to the particle-hole 
symmetry of graphene, the RKKY interaction induces 
ferromagnetic correlation between magnetic impurities 
on the same sublattice, while anti-ferromagnetic corre- 
lation between the ones on different sublattices. The de- 
pendence of the interaction on the distance R between 
two local magnetic moments, at the Dirac point, is found 
to be R~ 3 , whereas it behaves as R~ 2 in conventional 
two-dimensional (2D) systems^. Such a fast decay rate 
denotes that the interaction is rather short-ranged. In 
doped graphene, on the other hand, the spacial depen- 



dence of the interaction is predicted to be similar to con- 
ventional 2D systems, but this still remains to be exper- 
imentally verified. 

Due to the fact that the RKKY interaction is orig- 
inated by the exchange coupling between the impurity 
moments and the spin of itinerant electrons in the bulk 
of the system, spin polarization of electrons is expected to 
influence directly this interactional. In particular, combi- 
nation of the spin-dependence with a Dirac-like spectrum 
can mediate a much richer collective behavior of magnetic 
adatomsSS. This has been explained for surface states of 
a three dimensional topological insulator, on which mag- 
netic impurities exhibit a frustrated RKKY interaction 
with two possible phases of ordered ferromagnetic phase 
and a disordered spin glass phased Graphene, in partic- 
ular, with imbalanced chemical potentials of spin-up and 
spin-down electrons, presents a unique spin chiral ma- 
terial in which the interplay between the spin polariza- 
tion, gapless spectrum, and the chiral nature of electrons 
have been shown to result intriguing phenomenal. In 
two-dimension graphene system, the polarization of the 
chemical potential can be tuned to be of order or even 
higher than the mean chemical potential, the condition 
which is not possible in ordinary conductors. The aim of 
the present study is to address the question of how this 
peculiarity can affect the collective coupling of magnetic 
impurities on the surface of graphene sheet with a finite 
polarization of the spin. 

In this work, we calculate the RKKY interaction me- 
diated by spin-polarized Dirac fermions in a monolayer 
graphene using the Green's function method. Our theory 
for the spin polarization dependence of RKKY interac- 
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tion is motivated not only by fundamental transport con- 
siderations, but also by application and potential future 
experiments in graphene spintrotic field. With a spin- 
polarization along the z-axis, we show that the RKKY in- 
teraction is anisotropic corresponding to a XXZ model in- 
teraction between the two magnetic moments when their 
spin orientations are fixed. Besides R~ 3 dependence of 
the interaction for undoped graphene, we show particu- 
larly that the interaction behaves like R~ 2 when the spin 
polarization is finite. In addition, a beating pattern for 
the interaction in the cases where impurities are located 
along certain directions is obtained near the resonance 
condition which is controlled by the chemical potential 
and the spin polarization. 

The paper is organized as the following. In Sec. [TT| we 
introduce the formalism that will be used in calculating 
the RKKY interaction from the lattice Green's function. 
In Sec. |III| we present our analytic and numeric results for 
the coupling strengths of the RKKY interaction in both 
undoped and doped graphene sheets. Sec. |IV] contains 
discussions and a brief summary of our main results. 

II. METHOD AND THEORY 

We consider a spin-polarized graphene system identi- 
fied by a spin dependent chemical potential fi s (s = ±1), 
implying a mean chemical potential fi — ^ s fi s /2 and 
the spin polarization jj, p = ^2 s s^, s /2. Such a spin- 
polarization can be injected, for instance, by ferromag- 
netic electrodes on top of the graphene sheet 12212U_ Intrin- 
sic ferromagnetic correlations are also predicted to exist 
in graphene sheets^ and nanoribbons with zigzag edges 
^ under certain conditions. 

The electronic structure of spin polarized graphene can 
be reasonably well described using a rather simple tight- 
binding Hamiltonian, leading to analytical solutions for 
their energy dispersion and related eigenstates. The non- 
interacting nearest neighbor tight-binding Hamiltonian 
for 7r-band electrons with spin s, is determined by 2 ^ 

#5= - *E(4a>+&Ia>) 

(ij) 

~ s Mp {"*■-"'■- + & i,A«) . (!) 

i 

where a,i tS (bi lS ) annihilates an electron with spin s on 
sublattice A(B) of unit cell i and t ~ 2.9 eV denotes the 
nearest neighbor hopping parameter^. The sum (i, j) in 



Eq. ([I]) runs over distinct nearest neighbors. 

The s-component of the noninteracting Hamiltonian in 
momentum space is written as 

where the form factor in general case is /(k) = 
~*Ej e i m wn ich dj's are nearest neighbor posi- 
tion vectors. In this work, we are interested in the 
low-energy behavior, in which /(k) = vpk^(k), where 
$(fe) = e i^/3+o k ) at t h e Tji rac po i nt g an d = 
_ e »w3-0fe) a t t ne another Dirac point, K' , the chiral an- 
gle is Ok = tan~ 1 (k x /k y ), vf = 3ta/2h ~ 10 6 m/s is the 
Fermi velocity with a ~ 1.42A being the Carbon-Carbon 
distance in honeycomb lattice. 

Our system incorporates two localized magnetic mo- 
ments whose interaction is mediated through a spin po- 
larized electron liquid. We assume that the graphene is 
spin polarized first, and then we add the magnetic mo- 
ments. The contact interaction between the spin of itin- 
erant electrons and two magnetic impurities with mag- 
netic moments Mj and M 2 , located respectively at R x 
and R 2 , is given by 

H int = \J2 M r<Ki), (3) 

where A is the coupling constant between conduction elec- 
trons and impurity, s(r) = \ J2i H r i ~~ r ) <T i is the spin 
density operator^ with and <Tj being the position and 
vector of spin operators of ith electron. 

The RKKY interaction which arises from the quan- 
tum effects is obtained by using a second order pertur- 
batioiPEM] which reads as (from now on we set H = 1) 

_\ 2 roo 

#rkky = 3m / rf £ Tr[(M 1 . < r)G Q , 3 (R 1 ,R 2 ;e) 

J-oo 

x (M 2 -cr)G^(R 2 ,Ri;e)]n( M ) , (4) 

where n{fi) denotes the Fermi-Dirac distribution func- 
tion, <r is the vector of Pauli matrices in the spin space, 
G Q| g (Ri, R 2 ; e) is a 2 x 2 matrix of the single particle 
retarded Green's function in spin space, and a and j3 re- 
fer to the sublattices where two impurities are placed and 
finally, the trace is taken over the spin degree of freedom. 

For spin unpolarized graphene, Eq. Q simplifies to 
^rkky = jX( R i. R 2) Mi • M 2 , where x(Ri,R 2 ) is 
the spin susceptibility of the itinerant electrons and de- 
termines the indirect interaction between two local mo- 
ments. 
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In order to calculate the interaction Hamiltonian of 
Eq. Q, the form of the electronic single particle Green's 
function, G s (Ri,R 2 ;e) =< Ri|(e-M0+ -Wg)" 1 ^ > 
is needed. To calculate the retarded Green's function in 
real space, its Fourier components in momentum space 
might be first obtained. Due to the fact that our 2D Dirac 
fermion system is noninteracting and thus the direction 
of spin remains unchanged, the retarded Green's function 
G a p are diagonal in the spin space 



Gi A (R,0, £ ) = (e* K - 



' K H 1 ^'^)g AA {e- S ^) , (5) 



and 



X9ab(s - s^ip) 



(6) 



and 



Moreover, ^ BB - ^ AA 

G BA (0,R, e) = exp(— in/3)(exp(— iK ■ R — iOr) — 
exp(-iK' ■ R + i6 R ))g AB (e - s^ p ). Here 9aa{£) = 
^eK^—ieR/v-p) and jab(e) = JeKi(—ieR/vf), where 
K (x) and Ki(x) are the modified Bessel functions of 
the second kind, Or is the angle of the position R with 
respect to the K' — K direction and 7 = — 2ir/(Vlv ¥ ), in 
which n is the area of the Brillouin zone. 

By inserting the retarded Green's functions given by 
Eqs. Q and Q in Eq. Q, and taking the trace over spin 
degree of matrices, the RKKY Hamiltonian simplifies to 



H 



a/3 

RKKY 



= — [Jf{M lx M 2x +M ly M 2y ) + jfM u M 2z ] , 

(7) 

which is the honored XXZ model. Here J" 13 = 
-C<5> aP lf /R 3 and jf = -C<f> a/3 I? P / R 3 , with G = 
2(2n) 2 /fl 2 v F , <f> AA = 1 + cos[(K - K') • R] and = 
-1 + cos[(K - K') • R + 26 R ]. The different components 
of I a/3 for impurities on the same sublattice read as 



I AA = 2Qm 



jAA _ o. 



—00 

VF 



dy X-Ko(—ix-)x+Ko(—ix + ) 



s=± 



(8) 



where j/f = [aR/vf, x± = y ± Iif, y — sR/vf, and 
hp = ^pR/vF- For impurities on different sublattices, 
one only needs to replace Kq{x) with K\(x) in the above 
equations. Note that J BB = J AA and J BA = J AB . 

We find analytic results for the z component of the 
RKKY exchange coupling strength I z , for both cases 



where magnetic moments are located on the same or dif- 
ferent sublattices. For the same sublattice case, we begin 
by splitting the integral in the second line of Eq. ^ into 
the conduction and valance band contributions, and find 



jAA 



(R) = 



J-oo dx x 2 K 2 {~ix)^ 
E s =±^ dx x 2 K$(-ix)] 



(9) 

where xf± = 2/f±^f- The first integral can be solved^^ 
easily and the result is 7r 2 /32. The contribution from the 
second line of Eq. ([9|, can be also obtained by replac- 
ing Gm[ifQ(— ix)] with —ir 2 sign(x)Jo(\x\)Y (\x\)/2, and 
using the following relation 



dx x 2 sign(x) Jq(\x\)Y (\x\) 



Ff| 
2 X /^ 



M(x F ) , (10) 



where M(x) = G[{{i}, {§}}, {{1, 1}, 1}}, x 2 } is 
Meijer G-functiorP51. As a result, the function I AA is 
given by 



lf A (R) 



\XFs 



■M(x Fs 



(11) 



To calculate the long-range behavior of the RKKY 
interaction, the asymptotic behavior of the Mei- 
jer G-function is needed. It is also easy to see 
that asymptotic behavior of M(x) at large x is^ 
[2cos(2x) + 8xsin(2a;) — ir] /(8^/nx). It should be no- 
ticed that the M(x) tends to its long-range asymptotic 
expression for x > 2. Therefore, I A (R) for the long- 
range regime is simplified as 

I AA {R > a) ^ [cos(2x Fs ) + Ax Fs sin(2a: Fs )] 

(12) 

On the other hand, for the case that the impuri- 
ties are located on two different sublattices we can fol- 
low the same procedure discussed above, while we use 
J^dxx^mlKfi-ix)} = 3vr 2 /32, and %m[K 2 (-ix)\ = 



ir 2 sign(x)J 1 (\x\)Y 1 (\x\)/2 to find 



dxx 2 sign{x)J\ (\x\)Yi(\x\) 



2^ 



M'{x F ) , (13) 



where M'{x) = G[{{\}, {§}}, {{1, 2}, {-§, 0}}, x 2 ]. Fi- 
nally the I AB (R) reads as 



s=± 



M'(xfs) 



(14) 



The asymptotic behavior of M'(x) at large x is [37r 
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10cos(2x) - 8cc sin(2a:)]/(8 v ^r|a:|). Therefore, the long- 
range behavior of I AB is obtained as 

I AB {R^>a) « V [5cos(2x Fs ) +4x F;i sin(2x Fs )] . 
16 ' 

s=± 

(15) 

It should be mentioned that we were not able to find 
simple analytic expressions for the in plane components 
of the exchange coupling, and in the next section we 
will present our numerical results of them. 



III. NUMERICAL RESULTS AND 
DISCUSSIONS 



In this section, we present our main results for the 
RKKY exchange coupling in the presence of a spin po- 
larization Dirac fermions along the z-axis by analyzing 
the above calculated integrals of I x and I z . We extend 
the previously studied ^ results for dependencies on the 
distance R and lattice direction Br to the case of \i v ^ 0, 
for two different regimes of undoped, {ji = 0) and doped, 
(/i ^ 0) graphene. 

For the I z component of interactions, we solve the two 



expressions in Eqs. (11) and (14), numerically and then 



compare the results with asymptotic results obtained 



from the analytical expressions given by Eqs. (12 ) and 



(15), respectively. Generally, the results obtained from 



the two approaches match quite good in most of the case 
specially at long distances. The distance dependence of 
I z for both AA and AB cases are illustrated in Fig. [I] 
for the undoped graphene. For an unpolarized graphene 
fip = 0, I AA = I AA is just a constant. At finite /x p , 
the integral I z has quite different behavior exhibiting an 
oscillatory behavior as a function of R, with a linearly 
growing amplitude and a period given by 2ir/hp, as can 



be obtain directly from Eq. (12). This behavior of I AA 
results in an oscillatory J AA with a decreasing ampli- 
tude as R~ 2 , which mimics the behavior of the RKKY 
coupling of an unpolarized doped graphene^. We can 
understand this analogy by noting that the polarization 
induces spin-dependent doping of up and down spin Dirac 
bands of an undoped sample by shifting their chemical 
potential from the Dirac point. A comparison between 
independence of the integral I AA and that of I AB for 
various values of fj, p in Fig. |TJb) , shows their difference at 
short distance while reaching each other as R increases. 
As R — » the coupling interactions tend to their values 
of the unpolarized case where J AB is three times larger 
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FIG. 1: (Color online) (a) The integral I AA as a function of 
the distance R when both impurities are located on the same 
sublattice for various values of the spin polarization, /i p in 
units of eV. The chemical potential is set to zero. Symbols 
refer to the analytical results of Eq. (|12[ ) which are compared 
to the numerical evaluation of Eq. (1111), plotted as lines. For 



Mp = 0, h 



= H 



is just a constant, (b) A comparison 



between the integral I AA and I AB as a function of the distance 
R for various values of the spin polarization fi p in units of eV. 
At finite /j, p , the integral I z has a quite different behavior, 
oscillating as a function of R, with a period given by 2n/hF 
and a linearly growing amplitude. A comparison between in- 
dependence of the integral I AA and that of I AB , shows their 
difference at short distance while reaching each other as R 
increases. 



than J AB , as is discussed in Ref. [IB] . 

From the numerical calculations of the integrals ap- 
pearing in Eq. ([8]), we can also obtain the behavior of 
I AA and I AB for the RKKY interaction coupling of the 
components of the magnetic moments which are perpen- 
dicular to the spin polarization axis. Fig.^a) shows I AA 
as a function of R for the undoped graphene at different 



5 



< 

< X 



(a) 



H=0.0 

p 

— o— 11=0.1 

p 

H=0.5 

-Y-^1.0 y 

y 

▼ 

▼ 



StJOOQOOC 



,0000000000000000000 



R/a 




At finite values of both the chemical potential and 
the spin polarization, more complicated behavior of the 
RKKY coupling can be occurred. In this case, the be- 
havior of I z is determined by a superposition of four sinu- 
soidal functions with two different periods of 2it R/(xp—) 
and 2itR/ {xf+) each occurring twice with amplitudes 
l,Xp_ and l,Xp+, respectively. As the result, we ob- 
serve that for certain values of \i and fj, p , oscillations of 
I z exhibit a beating pattern with two characteristic pe- 
riods. Fig. |3ja) shows this beating behavior of integral, 
I A as a function of the impurities distance along arm- 
chair direction (where R = 3na with n being an integer 
number) for fi = 1.2 eV and fi p = 1.0 eV. Fig.pTb) shows 
the similar behavior of integral I AA as a distance along 
zigzag direction (where R = V3ma for an integer m) for 
/i = 1.2\/3 eV and /i p = V3 eV. We have obtained a sim- 
ilar beating pattern for oscillations of I AB , which also 
occurs for a certain values of \i and /i p . 

The behavior of the perpendicular components I x for 
fj, is also different from their linear behavior of the 
undoped case, as it is shown in Fig. Ha) for the fixed 
value of /j, = 0.5 eV and different values of /i p . In this 
case I X {R) exhibits oscillations with a linearly increas- 
ing amplitude whose slope is proportional to |/x — /i p |. 
Fig. Qb) is the same as Fig. |I]ja), but this time n P is 
fixed and \x changes. 



FIG. 2: (Color online) (a) The integral I AA versus distance 
between two impurities, when both impurities are located on 
the same sublattice for undoped graphene and for several val- 
ues of \i v in units of eV. Finite fi p , produces a linear increase 
of I AA with a slope proportional to /i p . (b) A comparison be- 
tween I x and I z for different configurations and for p p = 0.5 
eV. 



values of the spin polarization, /i p . For fj, p — 0, this 
function is a constant resulting a J AA which decays as 
R~ 3 . A finite difference fj, p between the chemical poten- 
tials of spin up and spin down carriers, produces a linear 
increase of I™ with a slope proportional to /i p . Thus, 
J AA decays as R~ 2 . Importantly, the sign of interaction 
J AA is always positive which shows that the coupling 
between the perpendicular components of the moments 
remains ferromagnetic-like for all Rs. To analyze the dif- 
ference between the two configurations of A A and AB, in 
Fig. [2^b) we have compared I AA and I AB , which shows 
that despite the difference at short distances, they tend 
to each other at larger distances. 



IV. SUMMARY AND CONCLUSIONS 

In conclusion, we have studied the influence of spin 
polarization on RKKY interaction in graphene. With 
a spin polarization along the z-axis, the induced inter- 
action between two magnetic impurities is found to be 
described by an anisotropic XXZ Hamiltonian with an 
exchange coupling depending on the distance R between 
the impurities and the doping level. For undoped but 
spin-polarized graphene, we have found that while the 
interaction between the x components of the moments 
remains constant with ferromagnetic sign, for the z com- 
ponents it oscillates with the distance R. In unpolarized 
spin case, the RKKY interaction induces ferromagnetic 
correlation between magnetic impurities on the same sub- 
lattice, while anti-ferromagnetic correlation between the 
ones on different sublattices^HHl. The dependence of the 
interaction on the distance R between two local magnetic 
moments, at the Dirac point, is found to be i? -3 , whereas 
it behaves as R~ 2 in doped graphene sheet. Besides R~ 3 
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J AB oc 7T 3 
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A* P = 


J AA oc - sin(2fc F i? + a)R-' 2 
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j£ A oc -~fipR~ 2 
J AA oc -sm(2k h R)iipR~ 2 
J AB oc ^pi?" 2 
J AS oc + sin(2fc h _R)/? -2 






J AA oc - sm{2k ¥ R)R-' z 
J AA oc -(sin(2a;p + ) + Csin(2:EF_)ir 2 
J AB oc sin(2fc F J?)i?" 2 
J Z AB oc (sin(2a;F+) + Csin(2a:F-)-R" 2 



TABLE I: A breakdown of the results on the scaling form of the RKKY interactions in monolayer graphitic systems. The 
results for vanishing spin polarization are reported from Ref. [16] while the ones for finite spin polarization denote to present 
work. It is worthwhile mentioning that the interaction Hamiltonian is modelled by the Heisenberg model for n P = and by 
the XX Z model for the case that n P ^ 0. We introduce a parameter fch = /-i p /vp and C — (/i — fi p )/(fj, + fj, p ) which are given 
by the chemical potential as well as the spin polarization. 



dependence of the interaction for undoped graphene, we 
show particularly that the interaction behaves like R~ 2 
when the spin polarization is finite. 

At finite value of both the chemical potential and 
the spin polarization, more complicated behavior of the 
RKKY coupling can be occurred. We have found that 
both components of the interaction oscillate with R. We 
have explored that for the chemical potentials fj, close to 
the polarization fi p , oscillations of the RKKY interaction 
exhibit a beating pattern when the impurities are located 
along zigzag or armchair directions. The two character- 
istic periods of the beating oscillations are determined by 



inverse of the difference and the sum of the chemical po- 
tential and the spin polarization. Since several works on 
RKKY interaction in 2D graphene systems are available, 
a proper comparison with those results seems to be in 
order (see Table I). 
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